##### for Extended Data Figure 10 a-b

rm(list=ls())
#加载包
library(Seurat)
library(dplyr)
library(tidyverse)
library(patchwork)
library(devtools)
library(harmony)
library(clustree)
library(fgsea)

##### data proccess #####
setwd("/Users/gin/fsdownload/spleen_sc/datadir")

dir_name <- list.files()

names(dir_name) = c('a20230321',"a20230328")  


scRNAlist <- list()
for(i in 1:length(dir_name)){
  counts <- Read10X(data.dir = dir_name[i],gene.column=1)
  scRNAlist[[i]] <- CreateSeuratObject(counts, min.cells = 50, min.features =200)
}

setwd("/Users/gin/fsdownload/spleen_sc/workdir")

for(i in 1:length(scRNAlist)){
  sc <- scRNAlist[[i]]
  sc[["mt_percent"]] <- PercentageFeatureSet(sc, pattern = "^MT-")
  HB_genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
  HB_m <- match(HB_genes, rownames(sc@assays$RNA))
  HB_genes <- rownames(sc@assays$RNA)[HB_m] 
  HB_genes <- HB_genes[!is.na(HB_genes)] 
  sc[["HB_percent"]] <- PercentageFeatureSet(sc, features=HB_genes) 
  scRNAlist[[i]] <- sc
  rm(sc)
}


scRNAlist <- lapply(X = scRNAlist, FUN = function(x){
  x <- subset(x, 
              subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & 
                mt_percent < 5 & 
                HB_percent < 3 & 
                nCount_RNA < quantile(nCount_RNA,0.97) & 
                nCount_RNA > 1000)})


scRNAlist <- merge(scRNAlist[[1]],y=c(scRNAlist[[2]]))


scRNAlist <- NormalizeData(scRNAlist) %>% 
  FindVariableFeatures(selection.method = "vst",nfeatures = 3000) %>% 
  ScaleData() %>% 
  RunPCA(npcs = 30, verbose = T)


scRNA_harmony <- RunHarmony(scRNAlist, group.by.vars = "orig.ident",project.dim = F)


scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:30)
scRNA_harmony <- FindClusters(object = scRNA_harmony,resolution = c(seq(0,1,.1)))
clustree(scRNA_harmony)
scRNA_harmony <- FindClusters(scRNA_harmony, resolution = 0.1)

scRNA_harmony <- RunTSNE(scRNA_harmony, reduction = "harmony", dims = 1:25)
scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:18)

DimPlot(scRNA_harmony,reduction = "tsne", pt.size = 1,label = T)
FeaturePlot(scRNA_harmony,features="GPC3",reduction="tsne",order=T,pt.size = 2,min.cutoff = 1)
FeaturePlot(scRNA_harmony,features="GREM1",reduction="tsne",order=T,pt.size = 2,min.cutoff = 1)
FeaturePlot(scRNA_harmony,features="MADCAM1",reduction="tsne",order=T,pt.size = 2,min.cutoff = 1)
FeaturePlot(scRNA_harmony,features="GPC3",reduction="tsne",order=T,pt.size = 2,min.cutoff = 1)


##### GSEA #####

library(presto)

genes = wilcoxauc(scRNA_harmony,"seurat_clusters")
gsea_genes<-genes %>% dplyr::filter(group=="1") %>% arrange(desc(auc)) %>% dplyr::select(feature,auc)
ranks = deframe(gsea_genes)


library(msigdbr)
m_df = msigdbr(species = "Homo sapiens", category = "C5")
fgsea_sets = m_df %>% split(x = .$gene_symbol, f = .$gs_name)

fgsea_res = fgsea(fgsea_sets, stats = ranks, nperm = 1000)

fgsea_tidy = fgsea_res %>% as_tibble() %>% arrange(desc(NES))

fgsea_tidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()

write.csv(fgsea_tidy[,1:7],"GSEA_C5.csv")

m_df = msigdbr(species = "Homo sapiens", category = "C2")
fgsea_sets = m_df %>% split(x = .$gene_symbol, f = .$gs_name)

fgsea_res = fgsea(fgsea_sets, stats = ranks, nperm = 1000)

fgsea_tidy = fgsea_res %>% as_tibble() %>% arrange(desc(NES))

fgsea_tidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% head()

write.csv(fgsea_tidy[,1:7],"GSEA_C2.csv")


##### membrane protein #####
memb = read.csv("protein_class_Predicted_membrane_proteins.tsv",sep="\t")
memb_list = unique(memb$Gene)

markers = FindMarkers(object = scRNA_harmony,ident.1 = 1,min.pct = 0.02,only.pos = T,logfc.threshold = 0.20)
markers$SYMBOL = rownames(markers)
markers_select = markers[markers$p_val_adj<0.01&markers$p_val<0.01&markers$avg_log2FC>0.20,]
markers_select = markers_select[markers_select$pct.2<0.01,]
write.csv(markers,"cluster8markers.csv")

marker_list = unique(markers_select$SYMBOL)

memb_p = intersect(marker_list,memb_list)
result = markers_select[memb_p,]

#FeaturePlot(scRNA_harmony,features="GREM1",reduction="tsne",order=T,pt.size = 1)
#VlnPlot(scRNA_harmony,features = result$SYMBOL[9],log=T)

write.csv(result,"cluster8memb2.csv")

fibro_memb = read.table("fibro_memb.txt",sep="\t",header=T)
fibro_memb_list = fibro_memb$Gene
tmp = intersect(fibro_memb_list,memb_p)

DimPlot(scRNA_harmony,reduction = "tsne",label = T,pt.size = 2)
FeaturePlot(scRNA_harmony,features="GPC3",reduction="tsne",order=T,pt.size = 2,min.cutoff = 1)



##### AddModuleScore #####
library(homologene)
library(AUCell)

mouse_gene = read.table("/Users/gin/fsdownload/mspleen_sc/allGrem1marker.csv",sep=",",header=T)
colnames(mouse_gene)[1] = "SYMBOL"
mouse_gene = mouse_gene[mouse_gene$p_val_adj<0.05,]
clustergene = mouse_gene$SYMBOL
clustergene = homologene(genes = clustergene,inTax = 10090,outTax = 9606)
clustergene = intersect(unique(clustergene[,2]),rownames(scRNA_harmony))
scRNA_harmony = AddModuleScore(scRNA_harmony,features = clustergene,name = "Grem1_score")
FeaturePlot(scRNA_harmony,features="Grem1_score1",reduction = "tsne",pt.size = 2,order=T,min.cutoff = 1)


##### singleR #####
library(SingleR)
library(celldex)

sce_for_SingleR <- GetAssayData(scRNA_harmony,slot = "counts")
clusters=scRNA_harmony@meta.data$seurat_clusters
mouseImmu <- ImmGenData()
pred.mouseImmu <- SingleR(test = sce_for_SingleR, ref = mouseImmu, labels = mouseImmu$label.main,
                          clusters = clusters, 
                          assay.type.test = "logcounts", assay.type.ref = "logcounts")


cellType=data.frame(ClusterID=levels(scRNA_harmony@meta.data$seurat_clusters),
                    mouseImmu=pred.mouseImmu$labels)
cellType



##### bubble plot #####
library(ggplot2)
bubble = read.table("selected GSEA.txt")
bubble = bubble[,c(1,2,5,7)]
bubble[,2] = -log10(bubble[,2])
#bubble[,4] = log10(bubble[,4])


labels=bubble[order(bubble$NES),"pathway"]
bubble$pathway = factor(bubble$pathway,levels=labels)

ggplot(bubble,aes(x=NES,y=pathway))+
  geom_point(aes(size=size,color=pval))+ 
  scale_colour_gradient(low="orange",high="red",)+
  scale_size_continuous(range = c(12,18)) + 
  xlim(c(1,5.5))+
  theme_bw()


